program mkxymaps
  use auty
  use suty
  use kuty
  implicit none
  type(ana_hdr) hdr
  interface 
     subroutine decs(e,xyz,etc)
       integer*8,intent(in)::e
       integer*4,intent(out)::xyz(3)
       integer*4,intent(out)::etc(4)
     end subroutine decs
  end interface      
  real,parameter::dx=2 ! mm
  real,parameter::dy=2
  real,parameter::dz=2 ! <<< PHONY
  !
  integer,parameter::nd=8   ! number of detectors
  integer::nx=512 ! 
  integer::ny=512
  !
  integer*4,target::xyz(3),etc(4)
  
  integer*4,allocatable::buf4(:)
  character*512 fin,fout,str
  integer*4 istat
  integer*8 n,ie,b1,b2,nb
  integer*8 fs
  integer*4,pointer::ix,iy
  integer jx,jy
  !
  logical lo_found
  integer nok,nlost,nnull,npair,nbtot
  integer*8 ev
  pointer(pev,ev)
  integer,allocatable::hist(:,:,:)
  integer ifile,nf
  integer istep
  integer pmtid
  integer narg,len

  !
  narg=command_argument_count()
  if(narg<4) then
     call usage
     call exit(1)
  end if

  call get_command_argument(1,str)
  nx=a2i(str)
  if(nx<=0) call claimer("*** Invalid nx: "//trim(str))

  call get_command_argument(2,str)
  ny=a2i(str)
  if(ny<=0) call claimer("*** Invalid ny: "//trim(str))

  call get_command_argument(3,fout,len)

  nf=narg-3

  ix=>xyz(1)
  iy=>xyz(2)

  ifile=0
  allocate(hist(nx,ny,nd))
  hist=0
  istep=0
  nbtot=0
  nok=0
  nlost=0
  nnull=0
  npair=0

  file: do ifile=4,narg
     call get_command_argument(ifile,fin)
     write(*,*) "Opening file: "//trim(fin)
     flush(6)
     open(unit=1,file=fin,status="old",access="stream",iostat=istat)
     if(istat/=0) call claimer("*** cannot open file: "//trim(fin))
     fs=fsize(fin)
     n=fs/4
     if(allocated(buf4)) deallocate(buf4)
     allocate(buf4(n),stat=istat)
     if(istat/=0) call claimer("*** cannot allocate buf")
     b1=1
     nb=1024*1024*1024/16
     
     list: do ie=1,n/nb+1
        b2=min(b1+nb,n)
        read(1,iostat=istat) buf4(b1:b2)    
        nbtot=nbtot+b2-b1+1
        b1=b2+1
        if(istat/=0) call claimer("*** error during read: "//trim(fin)//", code="//trim(i2a(istat)))
     end do list
     close(1)
     
     lo_found=.false.
     event: do ie=1,n
        if(buf4(ie)==0) then
           nnull=nnull+1
           cycle
        end if
        if(lo_found) then
           if(get_flg(buf4(ie))==0) then
              ! pair found
              npair=npair+1
              lo_found=.false.
              pev=loc(buf4(ie-1))
              call decs(ev,xyz,etc)
              if(etc(2)+1>nd) cycle
              pmtid=etc(2)+1
              jy=iy+1
              jx=ix+1
              if(jx>0.and.jx<=nx.and.jy>0.and.jy<=ny)&
                   hist(jx,jy,pmtid)=hist(jx,jy,pmtid)+1
              nok=nok+1
           else
              nlost=nlost+1
              lo_found=.true.
           end if
        else        
           if(get_flg(buf4(ie))==1) then
              lo_found=.true.
           else
              nlost=nlost+1
           end if
        end if
     end do event
  end do file

  call print_stat
  call write_hist(fout)

contains

  subroutine print_stat
    integer ip
    write(*,*) "Number of files read:      ",nf
    write(*,*) "Number of bytes read:      ",nb*4
    write(*,*) "Number of total events:    ",npair
    write(*,*) "Number of accept events:   ",nok
    write(*,*) "Number of lost 4-bytes:    ",nlost    
    do ip=1,nd
       write(*,*) "Block #"//trim(i2a(ip)),sum(hist(:,:,ip))
    end do    
  end subroutine print_stat

  integer function get_flg(e)
    integer*4,intent(in)::e
    get_flg=ishft(e,-30)
  end function get_flg
  
  subroutine write_header(fhdr)
    character*(*),intent(in)::fhdr
    call msets(loc(hdr),0,348)   
    hdr%hk%data_type="dsr"
    hdr%hk%db_name="hoge.dat"
    hdr%hk%sizeof_hdr=348
    hdr%hk%extents=16384
    hdr%hk%regular="r"
    hdr%id%dim(1)=4
    hdr%id%dim(2)=nx
    hdr%id%dim(3)=ny
    hdr%id%dim(4)=nd
    hdr%id%dim(5)=1
    hdr%id%pixdim(1)=4
    hdr%id%pixdim(2)=dx
    hdr%id%pixdim(3)=dy
    hdr%id%pixdim(4)=dz
    hdr%id%vox_units="mm"
    hdr%id%cal_units="count"
    hdr%id%datatype=DT_FLOAT
    hdr%id%bitpix=32
    hdr%id%cal_max=1
    hdr%id%glmax=1
    hdr%dh%generated="mkxymaps"
    hdr%dh%scannum="Unknown"
    hdr%dh%patient_id="Unknown"
    hdr%dh%descrip="2D-MAP"
    open(unit=3,file=fhdr,status="replace",access="stream",iostat=istat)
    if(istat/=0) call claimer("*** Cannot open file: "//trim(fhdr))
    write(3,iostat=istat) hdr
    if(istat/=0) call claimer("*** Error during write: "//trim(fhdr))
    close(3)
  end subroutine write_header
  
  subroutine write_hist(f)
    character*(*),intent(in)::f
    integer flen
    real*4 buf(nx,ny,nd)
    flen=len_trim(f) ! <<<    
    call write_header(f(1:flen-3)//"hdr")
    open(unit=3,file=f,iostat=istat,status="replace",access="stream")
    if(istat/=0) call claimer("*** cannot open file: "//trim(f))
    buf=hist
    buf=buf/maxval(buf)
    write(3,iostat=istat) buf
    close(3)
    if(istat/=0) call claimer("*** Error during write: "//trim(f))
  end subroutine write_hist
  
  subroutine claimer(msg,code)
    character*(*),intent(in),optional::msg
    integer,intent(in),optional::code
    integer rc
    if(present(msg)) write(0,"(x,a)") msg
    rc=-1
    if(present(code)) rc=code
    call exit(rc)
  end subroutine claimer

  subroutine usage
    write(*,100)
100 format(x, "Usage: mkxymaps nx ny output infiles...")
  end subroutine usage

end program mkxymaps
